% 这个脚本简易演示了弹簧振子模型，可以设置固有频率，阻力与周期性的外驱动力。
% 当存在阻力时，若欠阻尼（阻尼系数较小），振动振幅将指数减小，周期略微改变；若过阻尼（阻尼系数较大），振子将不振动
% 当存在阻力与外驱动力时，稳态后振子运动周期将与外驱动力相同；当外驱动力频率接近固有频率时振幅将增大
% CC BY 4

clc
clear
dt = 0.01;

k = 1; %弹簧弹性系数
m = 1; %振子质量
gamma = 0.5; %阻尼系数, 0 为无阻尼

Text = 6.28; %外驱动力周期
Fext_0 = 2; %外驱动力振幅, 0 为无外驱动力

w = sqrt(k/m);
T0 = 2*pi/w %固有周期，即无阻力也无外力时，振子的振动周期
wext = 2*pi/Text;

r(1) = 2; %初始位置
r(2) = 2;
t(1)=0;
t(2)=dt;

imgind = 1;

figure();
for i = 3:(10*floor(Text/dt))
  v = (r(i-1) -r(i-2))/(dt); %当前时刻振子速度
  f = - gamma*v; %阻力
  F = - k*r(i-1); %回复力（弹力）
  Fext = Fext_0*cos(wext*t(i-1)); %外驱动力

  r(i) = 2*r(i-1) - r(i-2) + (F+f+Fext)/m*(dt)^2;
  t(i) = (i-1)*dt;

  if mod(imgind,10) == 1 %绘制振子的运动动图
    clf;
    hold on
    axis equal
    axis ([-5 5 -5 5])

    y = r(i);

    scatter(0,0);
    scatter(0,y);

    quiver(0,y,0,0.2*f,'linewidth',0.8);
    quiver(0,y,0,0.2*Fext,'linewidth',0.8);
    quiver(0,y,0,0.2*F,'linewidth',0.8);
    drawnow;
    pause(0.01)
  endif
  imgind++;
end

figure() %绘制位移-时间图
hold on
plot(t, r);

xlabel('t')
ylabel('x')

